#include "defaults.h"
#include "nurbs.h"
#include "utils.h"
#include "triangle.h"
#include "triangleenergies.h"
#include "math/numerics.h"

namespace Optimizer {
static nfloat AbsDer(const nfloat& x) {
	return x < 0.0 ? -1.0 : 1.0;
}

#include "Mathematica.h"

nfloat BendingEnergyEdgeEvaluator::f(const Attribute& attribute) const {
	const EdgeAttribute& ea = static_cast<const EdgeAttribute&>(attribute);
	const boost::array<nfloat, 4>& x = ea.x;
	const boost::array<nfloat, 4>& y = ea.y;
	const boost::array<nfloat, 4>& z = ea.z;

	nfloat energy = 0.0;
	#include "CREnergy.inl"
	return energy;
}

Vector BendingEnergyEdgeEvaluator::Gradient(const Attribute& attribute) const {
	const EdgeAttribute& ea = static_cast<const EdgeAttribute&>(attribute);
	const boost::array<nfloat, 4>& x = ea.x;
	const boost::array<nfloat, 4>& y = ea.y;
	const boost::array<nfloat, 4>& z = ea.z;

	Vector grad(13);
	#include "CRGrad.inl"

	Vector grad2(12);
	for(size_t i = 0; i < 12; i++) grad2(i) = grad(i + 1);

	return grad2;
}

Matrix BendingEnergyEdgeEvaluator::Hessian(const Attribute& attribute) const {
	const EdgeAttribute& ea = static_cast<const EdgeAttribute&>(attribute);
	const boost::array<nfloat, 4>& x = ea.x;
	const boost::array<nfloat, 4>& y = ea.y;
	const boost::array<nfloat, 4>& z = ea.z;

	Matrix hess(13, 13);
	#include "CRHess.inl"

	Matrix hess2(12, 12);
	for(size_t i = 0; i < 12; i++)
		for(size_t j = 0; j < 12; j++) hess2(i, j) = hess(i + 1, j + 1);

	return hess2;
}

nfloat OrthogonalityEnergyTriangleEvaluator::f(const Attribute& attribute) const {
	const TriangleAttribute& ta = static_cast<const TriangleAttribute&>(attribute);
	STriangle* t = ta.t;
	
// 	SEdge* e = 0;
// 	if(     (t->edge[0]->neighbor == 0 && t->edge[1]->neighbor != 0 && t->edge[2]->neighbor != 0) ||
// 			(t->edge[0]->neighbor != 0 && t->edge[1]->neighbor == 0 && t->edge[2]->neighbor != 0) ||
// 			(t->edge[0]->neighbor != 0 && t->edge[1]->neighbor != 0 && t->edge[2]->neighbor == 0)) {
//  		for(size_t i = 0; i < 3; i++) {
// 			if(t->edge[i]->neighbor == 0)
// 				e = t->edge[i];
// 		}
// 	}
// 	if(e)
// 		return 0.0;
	
	const Vec3n& a = t->vertex[0]->Get<Vec3n>("position2");
	const Vec3n& b = t->vertex[1]->Get<Vec3n>("position2");
	const Vec3n& c = t->vertex[2]->Get<Vec3n>("position2");
	
	Matrix beta(3, 3);
	beta(0, 0) = a[0];
	beta(0, 1) = a[1];
	beta(0, 2) = 1.0;

	beta(1, 0) = b[0];
	beta(1, 1) = b[1];
	beta(1, 2) = 1.0;

	beta(2, 0) = c[0];
	beta(2, 1) = c[1];
	beta(2, 2) = 1.0;

	math::Invert(beta);
    
	Vec2n u, v;
	for(size_t i = 0; i < 3; i++) {
		const Vec2n gradi = math::V2n(beta(0, i), beta(1, i));
		const Vec3n param = t->vertex[i]->Get<Vec3n>(ta.strAttribute);
		u += param[0] * gradi;
		v += param[1] * gradi;
	}
	
	const nfloat lu = length(u);
	const nfloat lv = length(v);
	const nfloat eps = 0.001;

	if(lu > eps && lv > eps)
		return std::abs(Angle(u, v)) + std::abs(lu - lv);
	return 0.0;
}

Vector OrthogonalityEnergyTriangleEvaluator::Gradient(const Attribute& attribute) const {
	const TriangleAttribute& ta = static_cast<const TriangleAttribute&>(attribute);
	STriangle* t = ta.t;
	
// 	SEdge* e = 0;
// 	if(     (t->edge[0]->neighbor == 0 && t->edge[1]->neighbor != 0 && t->edge[2]->neighbor != 0) ||
// 			(t->edge[0]->neighbor != 0 && t->edge[1]->neighbor == 0 && t->edge[2]->neighbor != 0) ||
// 			(t->edge[0]->neighbor != 0 && t->edge[1]->neighbor != 0 && t->edge[2]->neighbor == 0)) {
//  		for(size_t i = 0; i < 3; i++) {
// 			if(t->edge[i]->neighbor == 0)
// 				e = t->edge[i];
// 		}
// 	}

	Vector g = ZeroVector(9);
// 	if(e)
// 		return g;

	
	const Vec3n& a = t->vertex[0]->Get<Vec3n>("position2");
	const Vec3n& b = t->vertex[1]->Get<Vec3n>("position2");
	const Vec3n& c = t->vertex[2]->Get<Vec3n>("position2");
	
	Matrix beta(3, 3);
	beta(0, 0) = a[0];
	beta(0, 1) = a[1];
	beta(0, 2) = 1.0;

	beta(1, 0) = b[0];
	beta(1, 1) = b[1];
	beta(1, 2) = 1.0;

	beta(2, 0) = c[0];
	beta(2, 1) = c[1];
	beta(2, 2) = 1.0;

	math::Invert(beta);
    
	Vec2n u, v;
	for(size_t i = 0; i < 3; i++) {
		const Vec2n gradi = math::V2n(beta(0, i), beta(1, i));
		const Vec3n param = t->vertex[i]->Get<Vec3n>(ta.strAttribute);
		u += param[0] * gradi;
		v += param[1] * gradi;
	}
	
	const Vec2n c1 = math::V2n(beta(0, 0), beta(1, 0));
	const Vec2n c2 = math::V2n(beta(0, 1), beta(1, 1));
	const Vec2n c3 = math::V2n(beta(0, 2), beta(1, 2));

	const nfloat cx1 = c1[0];
	const nfloat cy1 = c1[1];
	const nfloat cx2 = c2[0];
	const nfloat cy2 = c2[1];
	const nfloat cx3 = c3[0];
	const nfloat cy3 = c3[1];
	
	const nfloat ux = u[0];
	const nfloat uy = u[1];
	const nfloat vx = v[0];
	const nfloat vy = v[1];

	const nfloat lu2 = math::dot(u, u);
	const nfloat lv2 = math::dot(v, v);

	const nfloat lu = std::sqrt(lu2);
	const nfloat lv = std::sqrt(lv2);

	const nfloat eps = 0.001;
	if(lu > eps) {
		g(0) = (lu2*lv*(cx1*ux+cy1*uy)*AbsDer(lu-lv)-(cy1*(ux*uy*vx+(-lu2+std::pow(uy,2))*vy)+cx1*(-(lu2*vx)+ux*(ux*vx+uy*vy)))*AbsDer((ux*vx+uy*vy)/(lu*lv)))/(std::pow(lu2,1.5)*lv);
		g(1) = (-(lu*lv2*(cx1*vx+cy1*vy)*AbsDer(lu-lv))+(cx1*(lv2*ux-vx*(ux*vx+uy*vy))+cy1*(lv2*uy-vy*(ux*vx+uy*vy)))*AbsDer((ux*vx+uy*vy)/(lu*lv)))/(lu*std::pow(lv2,1.5));
		g(2) = 0.0;
		g(3) = (lu2*lv*(cx2*ux+cy2*uy)*AbsDer(lu-lv)-(cy2*(ux*uy*vx+(-lu2+std::pow(uy,2))*vy)+cx2*(-(lu2*vx)+ux*(ux*vx+uy*vy)))*AbsDer((ux*vx+uy*vy)/(lu*lv)))/(std::pow(lu2,1.5)*lv);
		g(4) = (-(lu*lv2*(cx2*vx+cy2*vy)*AbsDer(lu-lv))+(cx2*(lv2*ux-vx*(ux*vx+uy*vy))+cy2*(lv2*uy-vy*(ux*vx+uy*vy)))*AbsDer((ux*vx+uy*vy)/(lu*lv)))/(lu*std::pow(lv2,1.5));
		g(5) = 0.0;
		g(6) = (lu2*lv*(cx3*ux+cy3*uy)*AbsDer(lu-lv)-(cy3*(ux*uy*vx+(-lu2+std::pow(uy,2))*vy)+cx3*(-(lu2*vx)+ux*(ux*vx+uy*vy)))*AbsDer((ux*vx+uy*vy)/(lu*lv)))/(std::pow(lu2,1.5)*lv);
		g(7) = (-(lu*lv2*(cx3*vx+cy3*vy)*AbsDer(lu-lv))+(cx3*(lv2*ux-vx*(ux*vx+uy*vy))+cy3*(lv2*uy-vy*(ux*vx+uy*vy)))*AbsDer((ux*vx+uy*vy)/(lu*lv)))/(lu*std::pow(lv2,1.5));
		g(8) = 0.0;
	}

	return g;
}

Matrix OrthogonalityEnergyTriangleEvaluator::Hessian(const Attribute& attribute) const {
	return Matrix(0, 0);
}

nfloat BoundaryOrthogonalityEnergyTriangleEvaluator::f(const Attribute& attribute) const {
	const TriangleAttribute& ta = static_cast<const TriangleAttribute&>(attribute);
	STriangle* t = ta.t;
	
	const Vec3n& a = t->vertex[0]->Get<Vec3n>("position2");
	const Vec3n& b = t->vertex[1]->Get<Vec3n>("position2");
	const Vec3n& c = t->vertex[2]->Get<Vec3n>("position2");
	
	Matrix beta(3, 3);
	beta(0, 0) = a[0];
	beta(0, 1) = a[1];
	beta(0, 2) = 1.0;

	beta(1, 0) = b[0];
	beta(1, 1) = b[1];
	beta(1, 2) = 1.0;

	beta(2, 0) = c[0];
	beta(2, 1) = c[1];
	beta(2, 2) = 1.0;

	math::Invert(beta);
    
	Vec2n u, v;
	for(size_t i = 0; i < 3; i++) {
		const Vec2n gradi = math::V2n(beta(0, i), beta(1, i));
		const Vec3n param = t->vertex[i]->Get<Vec3n>(ta.strAttribute);
		u += param[0] * gradi;
		v += param[1] * gradi;
	}

	SEdge* e = 0;
	if(     (t->edge[0]->neighbor == 0 && t->edge[1]->neighbor != 0 && t->edge[2]->neighbor != 0) ||
			(t->edge[0]->neighbor != 0 && t->edge[1]->neighbor == 0 && t->edge[2]->neighbor != 0) ||
			(t->edge[0]->neighbor != 0 && t->edge[1]->neighbor != 0 && t->edge[2]->neighbor == 0)) {
 		for(size_t i = 0; i < 3; i++) {
			if(t->edge[i]->neighbor == 0)
				e = t->edge[i];
		}
	}

	if(e) {
		const Vec3n& a = e->a->Get<Vec3n>(ta.strAttribute);
		const Vec3n& b = e->b->Get<Vec3n>(ta.strAttribute);

		size_t type = 0;
		if(std::abs(a[0] - b[0]) < 0.00000001)
			type = 1;
		else
			type = 2;

		if(type == 0)
			std::cout << "isoparametric boundary edge" << std::endl;

		const Vec3n n = e->Normal("position2");	
		const Vec2n normal = math::V2n(n[0], n[1]);
		if(type == 1) {
			const nfloat Ev = std::abs(Angle(v, normal));
			return Ev;
		} else if(type == 2) {
			const nfloat Eu = std::abs(Angle(u, normal));
			return Eu;
		}
		
		const Vec3n center = e->Center("position2");
	}
	
	return 0.0;
}

Vector BoundaryOrthogonalityEnergyTriangleEvaluator::Gradient(const Attribute& attribute) const {
	const TriangleAttribute& ta = static_cast<const TriangleAttribute&>(attribute);
	STriangle* t = ta.t;
	
	const Vec3n& a = t->vertex[0]->Get<Vec3n>("position2");
	const Vec3n& b = t->vertex[1]->Get<Vec3n>("position2");
	const Vec3n& c = t->vertex[2]->Get<Vec3n>("position2");
	
	Matrix beta(3, 3);
	beta(0, 0) = a[0];
	beta(0, 1) = a[1];
	beta(0, 2) = 1.0;

	beta(1, 0) = b[0];
	beta(1, 1) = b[1];
	beta(1, 2) = 1.0;

	beta(2, 0) = c[0];
	beta(2, 1) = c[1];
	beta(2, 2) = 1.0;

	math::Invert(beta);
    
	Vec2n u, v;
	for(size_t i = 0; i < 3; i++) {
		const Vec2n gradi = math::V2n(beta(0, i), beta(1, i));
		const Vec3n param = t->vertex[i]->Get<Vec3n>(ta.strAttribute);
		u += param[0] * gradi;
		v += param[1] * gradi;
	}
	
	const Vec2n c1 = math::V2n(beta(0, 0), beta(1, 0));
	const Vec2n c2 = math::V2n(beta(0, 1), beta(1, 1));
	const Vec2n c3 = math::V2n(beta(0, 2), beta(1, 2));

	const nfloat cx1 = c1[0];
	const nfloat cy1 = c1[1];
	const nfloat cx2 = c2[0];
	const nfloat cy2 = c2[1];
	const nfloat cx3 = c3[0];
	const nfloat cy3 = c3[1];
	
	SEdge* e = 0;
	if(     (t->edge[0]->neighbor == 0 && t->edge[1]->neighbor != 0 && t->edge[2]->neighbor != 0) ||
			(t->edge[0]->neighbor != 0 && t->edge[1]->neighbor == 0 && t->edge[2]->neighbor != 0) ||
			(t->edge[0]->neighbor != 0 && t->edge[1]->neighbor != 0 && t->edge[2]->neighbor == 0)) {
 		for(size_t i = 0; i < 3; i++) {
			if(t->edge[i]->neighbor == 0)
				e = t->edge[i];
		}
	}

	Vector g = ZeroVector(9);
	if(e) {
		const Vec3n& a = e->a->Get<Vec3n>(ta.strAttribute);
		const Vec3n& b = e->b->Get<Vec3n>(ta.strAttribute);

		size_t type = 0;
		if(std::abs(a[0] - b[0]) < 0.00000001)
			type = 1;
		else
			type = 2;

		if(type == 0)
			std::cout << "isoparametric boundary edge" << std::endl;

		const Vec3n n = e->Normal("position2");	
		const Vec2n normal = math::V2n(n[0], n[1]);
		nfloat wx, wy, lw2;
		if(type == 1) {
			wx = v[0];
			wy = v[1];
			lw2 = dot(v, v);
		} else if(type == 2) {
			wx = u[0];
			wy = u[1];
			lw2 = dot(v, v);
		}
		const nfloat lw = std::sqrt(lw);
	
		const nfloat nx = normal[0];
		const nfloat ny = normal[1];
		const nfloat ln2 = math::dot(normal, normal);
		const nfloat ln = std::sqrt(ln2);

		if(type == 1) {
			g(0) = 0.0;
			g(1) = ((cx1*(lw2*nx-wx*(nx*wx+ny*wy))+cy1*(lw2*ny-wy*(nx*wx+ny*wy)))*AbsDer((nx*wx+ny*wy)/(ln*lw)))/(ln*Power(lw2,1.5));
			g(2) = 0.0;
			g(3) = 0.0;
			g(4) = ((cx2*(lw2*nx-wx*(nx*wx+ny*wy))+cy2*(lw2*ny-wy*(nx*wx+ny*wy)))*AbsDer((nx*wx+ny*wy)/(ln*lw)))/(ln*Power(lw2,1.5));
			g(5) = 0.0;
			g(6) = 0.0;
			g(7) = ((cx3*(lw2*nx-wx*(nx*wx+ny*wy))+cy3*(lw2*ny-wy*(nx*wx+ny*wy)))*AbsDer((nx*wx+ny*wy)/(ln*lw)))/(ln*Power(lw2,1.5));
			g(8) = 0.0;
		} else if(type == 2) {
			g(0) = ((cx1*(lw2*nx-wx*(nx*wx+ny*wy))+cy1*(lw2*ny-wy*(nx*wx+ny*wy)))*AbsDer((nx*wx+ny*wy)/(ln*lw)))/(ln*Power(lw2,1.5));
			g(1) = 0.0;
			g(2) = 0.0;
			g(3) = ((cx2*(lw2*nx-wx*(nx*wx+ny*wy))+cy2*(lw2*ny-wy*(nx*wx+ny*wy)))*AbsDer((nx*wx+ny*wy)/(ln*lw)))/(ln*Power(lw2,1.5));
			g(4) = 0.0;
			g(5) = 0.0;
			g(6) = ((cx3*(lw2*nx-wx*(nx*wx+ny*wy))+cy3*(lw2*ny-wy*(nx*wx+ny*wy)))*AbsDer((nx*wx+ny*wy)/(ln*lw)))/(ln*Power(lw2,1.5));
			g(7) = 0.0;
			g(8) = 0.0;
		}
	}

	return g;
}

Matrix BoundaryOrthogonalityEnergyTriangleEvaluator::Hessian(const Attribute& attribute) const {
	return Matrix(0, 0);
}

nfloat IsometricEnergyTriangleEvaluator::f(const Attribute& attribute) const {
	const TriangleAttribute& ta = static_cast<const TriangleAttribute&>(attribute);
	STriangle* t = ta.t;

	const Vec3n& v1 = t->vertex[0]->Get<Vec3n>(ta.strAttribute);
	const Vec3n& v2 = t->vertex[1]->Get<Vec3n>(ta.strAttribute);
	const Vec3n& v3 = t->vertex[2]->Get<Vec3n>(ta.strAttribute);
	
	const Vec3n& b = t->Get<Vec3n>("beta");
	const Vec3n& w = t->Get<Vec3n>("weight");

	Vec3n a1 = v2 - v1;
	Vec3n b1 = v3 - v1;

	Vec3n a2 = v1 - v2;
	Vec3n b2 = v3 - v2;

	Vec3n a3 = v2 - v3;
	Vec3n b3 = v1 - v3;

	normalize(a1);
	normalize(b1);
	normalize(a2);
	normalize(b2);
	normalize(a3);
	normalize(b3);

	const Vec3n a = math::V3n(std::acos(dot(a1, b1)), std::acos(dot(a2, b2)), std::acos(dot(a3, b3)));

	nfloat E = 0.0;
	for(size_t i = 0; i < 3; i++)
		E += (1.0 / w[i]) * ((a[i] - b[i]) * (a[i] - b[i])); 
	
	return E;
}

Vector IsometricEnergyTriangleEvaluator::Gradient(const Attribute& attribute) const {
	const TriangleAttribute& ta = static_cast<const TriangleAttribute&>(attribute);
	STriangle* t = ta.t;

	const Vec3n& p1 = t->vertex[0]->Get<Vec3n>(ta.strAttribute);
	const Vec3n& p2 = t->vertex[1]->Get<Vec3n>(ta.strAttribute);
	const Vec3n& p3 = t->vertex[2]->Get<Vec3n>(ta.strAttribute);
	
	const Vec3n& b = t->Get<Vec3n>("beta");
	const Vec3n& w = t->Get<Vec3n>("weight");

	const nfloat x1 = p1[0];
	const nfloat x2 = p2[0];
	const nfloat x3 = p3[0];
	const nfloat y1 = p1[1];
	const nfloat y2 = p2[1];
	const nfloat y3 = p3[1];
	const nfloat z1 = p1[2];
	const nfloat z2 = p2[2];
	const nfloat z3 = p3[2];
	const nfloat w1 = w[0];
	const nfloat w2 = w[1];
	const nfloat w3 = w[2];
	const nfloat b1 = b[0];
	const nfloat b2 = b[1];
	const nfloat b3 = b[2];

	const nfloat s1=std::sqrt(std::pow(std::abs(x1 - x2),2) + std::pow(std::abs(y1 - y2),2) + std::pow(std::abs(z1 - z2),2));
	const nfloat s2=std::sqrt(std::pow(std::abs(-x2 + x3),2) + std::pow(std::abs(-y2 + y3),2) + std::pow(std::abs(-z2 + z3),2));
	const nfloat s3=std::pow(std::pow(std::abs(x1 - x2),2) + std::pow(std::abs(y1 - y2),2) + std::pow(std::abs(z1 - z2),2),1.5);
	const nfloat s4=std::sqrt(std::pow(std::abs(x1 - x3),2) + std::pow(std::abs(y1 - y3),2) + std::pow(std::abs(z1 - z3),2));
	const nfloat s5=std::sqrt(std::pow(std::abs(x2 - x3),2) + std::pow(std::abs(y2 - y3),2) + std::pow(std::abs(z2 - z3),2));
	const nfloat s6=std::pow(std::pow(std::abs(x1 - x3),2) + std::pow(std::abs(y1 - y3),2) + std::pow(std::abs(z1 - z3),2),1.5);
	const nfloat s7=std::sqrt(std::pow(std::abs(-x1 + x2),2) + std::pow(std::abs(-y1 + y2),2) + std::pow(std::abs(-z1 + z2),2));
	const nfloat s8=std::sqrt(std::pow(std::abs(-x1 + x3),2) + std::pow(std::abs(-y1 + y3),2) + std::pow(std::abs(-z1 + z3),2));
	const nfloat s9=std::pow(std::pow(std::abs(-x1 + x2),2) + std::pow(std::abs(-y1 + y2),2) + std::pow(std::abs(-z1 + z2),2),1.5);
	const nfloat s10=std::pow(std::pow(std::abs(x2 - x3),2) + std::pow(std::abs(y2 - y3),2) + std::pow(std::abs(z2 - z3),2),1.5);
	const nfloat s11=std::pow(std::pow(std::abs(-x1 + x3),2) + std::pow(std::abs(-y1 + y3),2) + std::pow(std::abs(-z1 + z3),2),1.5);
	const nfloat s12=std::pow(std::pow(std::abs(-x2 + x3),2) + std::pow(std::abs(-y2 + y3),2) + std::pow(std::abs(-z2 + z3),2),1.5);

	Vector g(9);
	g(0)=(-2*(-b2+std::acos(((x1-x2)*(-x2+x3))/(s1*s2)+((y1-y2)*(-y2+y3))/(s1*s2)+((z1-z2)*(-z2+z3))/(s1*s2)))*((-x2+x3)/(s1*s2)-((x1-x2)*(-x2+x3)*std::abs(x1-x2)*AbsDer(x1-x2))/(s3*s2)-((y1-y2)*(-y2+y3)*std::abs(x1-x2)*AbsDer(x1-x2))/(s3*s2)-((z1-z2)*(-z2+z3)*std::abs(x1-x2)*AbsDer(x1-x2))/(s3*s2)))/(w2*std::sqrt(1-std::pow(((x1-x2)*(-x2+x3))/(s1*s2)+((y1-y2)*(-y2+y3))/(s1*s2)+((z1-z2)*(-z2+z3))/(s1*s2),2)))-(2*(-b3+std::acos(((x1-x3)*(x2-x3))/(s4*s5)+((y1-y3)*(y2-y3))/(s4*s5)+((z1-z3)*(z2-z3))/(s4*s5)))*((x2-x3)/(s4*s5)-((x1-x3)*(x2-x3)*std::abs(x1-x3)*AbsDer(x1-x3))/(s6*s5)-((y1-y3)*(y2-y3)*std::abs(x1-x3)*AbsDer(x1-x3))/(s6*s5)-((z1-z3)*(z2-z3)*std::abs(x1-x3)*AbsDer(x1-x3))/(s6*s5)))/(w3*std::sqrt(1-std::pow(((x1-x3)*(x2-x3))/(s4*s5)+((y1-y3)*(y2-y3))/(s4*s5)+((z1-z3)*(z2-z3))/(s4*s5),2)))-(2*(-b1+std::acos(((-x1+x2)*(-x1+x3))/(s7*s8)+((-y1+y2)*(-y1+y3))/(s7*s8)+((-z1+z2)*(-z1+z3))/(s7*s8)))*(-((-x1+x2)/(s7*s8))-(-x1+x3)/(s7*s8)+((-x1+x2)*(-x1+x3)*std::abs(-x1+x2)*AbsDer(-x1+x2))/(s9*s8)+((-y1+y2)*(-y1+y3)*std::abs(-x1+x2)*AbsDer(-x1+x2))/(s9*s8)+((-z1+z2)*(-z1+z3)*std::abs(-x1+x2)*AbsDer(-x1+x2))/(s9*s8)+((-x1+x2)*(-x1+x3)*std::abs(-x1+x3)*AbsDer(-x1+x3))/(s7*s11)+((-y1+y2)*(-y1+y3)*std::abs(-x1+x3)*AbsDer(-x1+x3))/(s7*s11)+((-z1+z2)*(-z1+z3)*std::abs(-x1+x3)*AbsDer(-x1+x3))/(s7*s11)))/(w1*std::sqrt(1-std::pow(((-x1+x2)*(-x1+x3))/(s7*s8)+((-y1+y2)*(-y1+y3))/(s7*s8)+((-z1+z2)*(-z1+z3))/(s7*s8),2)));
	g(1)=(-2*(-b2+std::acos(((x1-x2)*(-x2+x3))/(s1*s2)+((y1-y2)*(-y2+y3))/(s1*s2)+((z1-z2)*(-z2+z3))/(s1*s2)))*((-y2+y3)/(s1*s2)-((x1-x2)*(-x2+x3)*std::abs(y1-y2)*AbsDer(y1-y2))/(s3*s2)-((y1-y2)*(-y2+y3)*std::abs(y1-y2)*AbsDer(y1-y2))/(s3*s2)-((z1-z2)*(-z2+z3)*std::abs(y1-y2)*AbsDer(y1-y2))/(s3*s2)))/(w2*std::sqrt(1-std::pow(((x1-x2)*(-x2+x3))/(s1*s2)+((y1-y2)*(-y2+y3))/(s1*s2)+((z1-z2)*(-z2+z3))/(s1*s2),2)))-(2*(-b3+std::acos(((x1-x3)*(x2-x3))/(s4*s5)+((y1-y3)*(y2-y3))/(s4*s5)+((z1-z3)*(z2-z3))/(s4*s5)))*((y2-y3)/(s4*s5)-((x1-x3)*(x2-x3)*std::abs(y1-y3)*AbsDer(y1-y3))/(s6*s5)-((y1-y3)*(y2-y3)*std::abs(y1-y3)*AbsDer(y1-y3))/(s6*s5)-((z1-z3)*(z2-z3)*std::abs(y1-y3)*AbsDer(y1-y3))/(s6*s5)))/(w3*std::sqrt(1-std::pow(((x1-x3)*(x2-x3))/(s4*s5)+((y1-y3)*(y2-y3))/(s4*s5)+((z1-z3)*(z2-z3))/(s4*s5),2)))-(2*(-b1+std::acos(((-x1+x2)*(-x1+x3))/(s7*s8)+((-y1+y2)*(-y1+y3))/(s7*s8)+((-z1+z2)*(-z1+z3))/(s7*s8)))*(-((-y1+y2)/(s7*s8))-(-y1+y3)/(s7*s8)+((-x1+x2)*(-x1+x3)*std::abs(-y1+y2)*AbsDer(-y1+y2))/(s9*s8)+((-y1+y2)*(-y1+y3)*std::abs(-y1+y2)*AbsDer(-y1+y2))/(s9*s8)+((-z1+z2)*(-z1+z3)*std::abs(-y1+y2)*AbsDer(-y1+y2))/(s9*s8)+((-x1+x2)*(-x1+x3)*std::abs(-y1+y3)*AbsDer(-y1+y3))/(s7*s11)+((-y1+y2)*(-y1+y3)*std::abs(-y1+y3)*AbsDer(-y1+y3))/(s7*s11)+((-z1+z2)*(-z1+z3)*std::abs(-y1+y3)*AbsDer(-y1+y3))/(s7*s11)))/(w1*std::sqrt(1-std::pow(((-x1+x2)*(-x1+x3))/(s7*s8)+((-y1+y2)*(-y1+y3))/(s7*s8)+((-z1+z2)*(-z1+z3))/(s7*s8),2)));
	g(2)=(-2*(-b2+std::acos(((x1-x2)*(-x2+x3))/(s1*s2)+((y1-y2)*(-y2+y3))/(s1*s2)+((z1-z2)*(-z2+z3))/(s1*s2)))*((-z2+z3)/(s1*s2)-((x1-x2)*(-x2+x3)*std::abs(z1-z2)*AbsDer(z1-z2))/(s3*s2)-((y1-y2)*(-y2+y3)*std::abs(z1-z2)*AbsDer(z1-z2))/(s3*s2)-((z1-z2)*(-z2+z3)*std::abs(z1-z2)*AbsDer(z1-z2))/(s3*s2)))/(w2*std::sqrt(1-std::pow(((x1-x2)*(-x2+x3))/(s1*s2)+((y1-y2)*(-y2+y3))/(s1*s2)+((z1-z2)*(-z2+z3))/(s1*s2),2)))-(2*(-b3+std::acos(((x1-x3)*(x2-x3))/(s4*s5)+((y1-y3)*(y2-y3))/(s4*s5)+((z1-z3)*(z2-z3))/(s4*s5)))*((z2-z3)/(s4*s5)-((x1-x3)*(x2-x3)*std::abs(z1-z3)*AbsDer(z1-z3))/(s6*s5)-((y1-y3)*(y2-y3)*std::abs(z1-z3)*AbsDer(z1-z3))/(s6*s5)-((z1-z3)*(z2-z3)*std::abs(z1-z3)*AbsDer(z1-z3))/(s6*s5)))/(w3*std::sqrt(1-std::pow(((x1-x3)*(x2-x3))/(s4*s5)+((y1-y3)*(y2-y3))/(s4*s5)+((z1-z3)*(z2-z3))/(s4*s5),2)))-(2*(-b1+std::acos(((-x1+x2)*(-x1+x3))/(s7*s8)+((-y1+y2)*(-y1+y3))/(s7*s8)+((-z1+z2)*(-z1+z3))/(s7*s8)))*(-((-z1+z2)/(s7*s8))-(-z1+z3)/(s7*s8)+((-x1+x2)*(-x1+x3)*std::abs(-z1+z2)*AbsDer(-z1+z2))/(s9*s8)+((-y1+y2)*(-y1+y3)*std::abs(-z1+z2)*AbsDer(-z1+z2))/(s9*s8)+((-z1+z2)*(-z1+z3)*std::abs(-z1+z2)*AbsDer(-z1+z2))/(s9*s8)+((-x1+x2)*(-x1+x3)*std::abs(-z1+z3)*AbsDer(-z1+z3))/(s7*s11)+((-y1+y2)*(-y1+y3)*std::abs(-z1+z3)*AbsDer(-z1+z3))/(s7*s11)+((-z1+z2)*(-z1+z3)*std::abs(-z1+z3)*AbsDer(-z1+z3))/(s7*s11)))/(w1*std::sqrt(1-std::pow(((-x1+x2)*(-x1+x3))/(s7*s8)+((-y1+y2)*(-y1+y3))/(s7*s8)+((-z1+z2)*(-z1+z3))/(s7*s8),2)));
	g(3)=(-2*(-b1+std::acos(((-x1+x2)*(-x1+x3))/(s7*s8)+((-y1+y2)*(-y1+y3))/(s7*s8)+((-z1+z2)*(-z1+z3))/(s7*s8)))*((-x1+x3)/(s7*s8)-((-x1+x2)*(-x1+x3)*std::abs(-x1+x2)*AbsDer(-x1+x2))/(s9*s8)-((-y1+y2)*(-y1+y3)*std::abs(-x1+x2)*AbsDer(-x1+x2))/(s9*s8)-((-z1+z2)*(-z1+z3)*std::abs(-x1+x2)*AbsDer(-x1+x2))/(s9*s8)))/(w1*std::sqrt(1-std::pow(((-x1+x2)*(-x1+x3))/(s7*s8)+((-y1+y2)*(-y1+y3))/(s7*s8)+((-z1+z2)*(-z1+z3))/(s7*s8),2)))-(2*(-b3+std::acos(((x1-x3)*(x2-x3))/(s4*s5)+((y1-y3)*(y2-y3))/(s4*s5)+((z1-z3)*(z2-z3))/(s4*s5)))*((x1-x3)/(s4*s5)-((x1-x3)*(x2-x3)*std::abs(x2-x3)*AbsDer(x2-x3))/(s4*s10)-((y1-y3)*(y2-y3)*std::abs(x2-x3)*AbsDer(x2-x3))/(s4*s10)-((z1-z3)*(z2-z3)*std::abs(x2-x3)*AbsDer(x2-x3))/(s4*s10)))/(w3*std::sqrt(1-std::pow(((x1-x3)*(x2-x3))/(s4*s5)+((y1-y3)*(y2-y3))/(s4*s5)+((z1-z3)*(z2-z3))/(s4*s5),2)))-(2*(-b2+std::acos(((x1-x2)*(-x2+x3))/(s1*s2)+((y1-y2)*(-y2+y3))/(s1*s2)+((z1-z2)*(-z2+z3))/(s1*s2)))*(-((x1-x2)/(s1*s2))-(-x2+x3)/(s1*s2)+((x1-x2)*(-x2+x3)*std::abs(x1-x2)*AbsDer(x1-x2))/(s3*s2)+((y1-y2)*(-y2+y3)*std::abs(x1-x2)*AbsDer(x1-x2))/(s3*s2)+((z1-z2)*(-z2+z3)*std::abs(x1-x2)*AbsDer(x1-x2))/(s3*s2)+((x1-x2)*(-x2+x3)*std::abs(-x2+x3)*AbsDer(-x2+x3))/(s1*s12)+((y1-y2)*(-y2+y3)*std::abs(-x2+x3)*AbsDer(-x2+x3))/(s1*s12)+((z1-z2)*(-z2+z3)*std::abs(-x2+x3)*AbsDer(-x2+x3))/(s1*s12)))/(w2*std::sqrt(1-std::pow(((x1-x2)*(-x2+x3))/(s1*s2)+((y1-y2)*(-y2+y3))/(s1*s2)+((z1-z2)*(-z2+z3))/(s1*s2),2)));
	g(4)=(-2*(-b1+std::acos(((-x1+x2)*(-x1+x3))/(s7*s8)+((-y1+y2)*(-y1+y3))/(s7*s8)+((-z1+z2)*(-z1+z3))/(s7*s8)))*((-y1+y3)/(s7*s8)-((-x1+x2)*(-x1+x3)*std::abs(-y1+y2)*AbsDer(-y1+y2))/(s9*s8)-((-y1+y2)*(-y1+y3)*std::abs(-y1+y2)*AbsDer(-y1+y2))/(s9*s8)-((-z1+z2)*(-z1+z3)*std::abs(-y1+y2)*AbsDer(-y1+y2))/(s9*s8)))/(w1*std::sqrt(1-std::pow(((-x1+x2)*(-x1+x3))/(s7*s8)+((-y1+y2)*(-y1+y3))/(s7*s8)+((-z1+z2)*(-z1+z3))/(s7*s8),2)))-(2*(-b3+std::acos(((x1-x3)*(x2-x3))/(s4*s5)+((y1-y3)*(y2-y3))/(s4*s5)+((z1-z3)*(z2-z3))/(s4*s5)))*((y1-y3)/(s4*s5)-((x1-x3)*(x2-x3)*std::abs(y2-y3)*AbsDer(y2-y3))/(s4*s10)-((y1-y3)*(y2-y3)*std::abs(y2-y3)*AbsDer(y2-y3))/(s4*s10)-((z1-z3)*(z2-z3)*std::abs(y2-y3)*AbsDer(y2-y3))/(s4*s10)))/(w3*std::sqrt(1-std::pow(((x1-x3)*(x2-x3))/(s4*s5)+((y1-y3)*(y2-y3))/(s4*s5)+((z1-z3)*(z2-z3))/(s4*s5),2)))-(2*(-b2+std::acos(((x1-x2)*(-x2+x3))/(s1*s2)+((y1-y2)*(-y2+y3))/(s1*s2)+((z1-z2)*(-z2+z3))/(s1*s2)))*(-((y1-y2)/(s1*s2))-(-y2+y3)/(s1*s2)+((x1-x2)*(-x2+x3)*std::abs(y1-y2)*AbsDer(y1-y2))/(s3*s2)+((y1-y2)*(-y2+y3)*std::abs(y1-y2)*AbsDer(y1-y2))/(s3*s2)+((z1-z2)*(-z2+z3)*std::abs(y1-y2)*AbsDer(y1-y2))/(s3*s2)+((x1-x2)*(-x2+x3)*std::abs(-y2+y3)*AbsDer(-y2+y3))/(s1*s12)+((y1-y2)*(-y2+y3)*std::abs(-y2+y3)*AbsDer(-y2+y3))/(s1*s12)+((z1-z2)*(-z2+z3)*std::abs(-y2+y3)*AbsDer(-y2+y3))/(s1*s12)))/(w2*std::sqrt(1-std::pow(((x1-x2)*(-x2+x3))/(s1*s2)+((y1-y2)*(-y2+y3))/(s1*s2)+((z1-z2)*(-z2+z3))/(s1*s2),2)));
	g(5)=(-2*(-b1+std::acos(((-x1+x2)*(-x1+x3))/(s7*s8)+((-y1+y2)*(-y1+y3))/(s7*s8)+((-z1+z2)*(-z1+z3))/(s7*s8)))*((-z1+z3)/(s7*s8)-((-x1+x2)*(-x1+x3)*std::abs(-z1+z2)*AbsDer(-z1+z2))/(s9*s8)-((-y1+y2)*(-y1+y3)*std::abs(-z1+z2)*AbsDer(-z1+z2))/(s9*s8)-((-z1+z2)*(-z1+z3)*std::abs(-z1+z2)*AbsDer(-z1+z2))/(s9*s8)))/(w1*std::sqrt(1-std::pow(((-x1+x2)*(-x1+x3))/(s7*s8)+((-y1+y2)*(-y1+y3))/(s7*s8)+((-z1+z2)*(-z1+z3))/(s7*s8),2)))-(2*(-b3+std::acos(((x1-x3)*(x2-x3))/(s4*s5)+((y1-y3)*(y2-y3))/(s4*s5)+((z1-z3)*(z2-z3))/(s4*s5)))*((z1-z3)/(s4*s5)-((x1-x3)*(x2-x3)*std::abs(z2-z3)*AbsDer(z2-z3))/(s4*s10)-((y1-y3)*(y2-y3)*std::abs(z2-z3)*AbsDer(z2-z3))/(s4*s10)-((z1-z3)*(z2-z3)*std::abs(z2-z3)*AbsDer(z2-z3))/(s4*s10)))/(w3*std::sqrt(1-std::pow(((x1-x3)*(x2-x3))/(s4*s5)+((y1-y3)*(y2-y3))/(s4*s5)+((z1-z3)*(z2-z3))/(s4*s5),2)))-(2*(-b2+std::acos(((x1-x2)*(-x2+x3))/(s1*s2)+((y1-y2)*(-y2+y3))/(s1*s2)+((z1-z2)*(-z2+z3))/(s1*s2)))*(-((z1-z2)/(s1*s2))-(-z2+z3)/(s1*s2)+((x1-x2)*(-x2+x3)*std::abs(z1-z2)*AbsDer(z1-z2))/(s3*s2)+((y1-y2)*(-y2+y3)*std::abs(z1-z2)*AbsDer(z1-z2))/(s3*s2)+((z1-z2)*(-z2+z3)*std::abs(z1-z2)*AbsDer(z1-z2))/(s3*s2)+((x1-x2)*(-x2+x3)*std::abs(-z2+z3)*AbsDer(-z2+z3))/(s1*s12)+((y1-y2)*(-y2+y3)*std::abs(-z2+z3)*AbsDer(-z2+z3))/(s1*s12)+((z1-z2)*(-z2+z3)*std::abs(-z2+z3)*AbsDer(-z2+z3))/(s1*s12)))/(w2*std::sqrt(1-std::pow(((x1-x2)*(-x2+x3))/(s1*s2)+((y1-y2)*(-y2+y3))/(s1*s2)+((z1-z2)*(-z2+z3))/(s1*s2),2)));
	g(6)=(-2*(-b3+std::acos(((x1-x3)*(x2-x3))/(s4*s5)+((y1-y3)*(y2-y3))/(s4*s5)+((z1-z3)*(z2-z3))/(s4*s5)))*(-((x1-x3)/(s4*s5))-(x2-x3)/(s4*s5)+((x1-x3)*(x2-x3)*std::abs(x1-x3)*AbsDer(x1-x3))/(s6*s5)+((y1-y3)*(y2-y3)*std::abs(x1-x3)*AbsDer(x1-x3))/(s6*s5)+((z1-z3)*(z2-z3)*std::abs(x1-x3)*AbsDer(x1-x3))/(s6*s5)+((x1-x3)*(x2-x3)*std::abs(x2-x3)*AbsDer(x2-x3))/(s4*s10)+((y1-y3)*(y2-y3)*std::abs(x2-x3)*AbsDer(x2-x3))/(s4*s10)+((z1-z3)*(z2-z3)*std::abs(x2-x3)*AbsDer(x2-x3))/(s4*s10)))/(w3*std::sqrt(1-std::pow(((x1-x3)*(x2-x3))/(s4*s5)+((y1-y3)*(y2-y3))/(s4*s5)+((z1-z3)*(z2-z3))/(s4*s5),2)))-(2*(-b1+std::acos(((-x1+x2)*(-x1+x3))/(s7*s8)+((-y1+y2)*(-y1+y3))/(s7*s8)+((-z1+z2)*(-z1+z3))/(s7*s8)))*((-x1+x2)/(s7*s8)-((-x1+x2)*(-x1+x3)*std::abs(-x1+x3)*AbsDer(-x1+x3))/(s7*s11)-((-y1+y2)*(-y1+y3)*std::abs(-x1+x3)*AbsDer(-x1+x3))/(s7*s11)-((-z1+z2)*(-z1+z3)*std::abs(-x1+x3)*AbsDer(-x1+x3))/(s7*s11)))/(w1*std::sqrt(1-std::pow(((-x1+x2)*(-x1+x3))/(s7*s8)+((-y1+y2)*(-y1+y3))/(s7*s8)+((-z1+z2)*(-z1+z3))/(s7*s8),2)))-(2*(-b2+std::acos(((x1-x2)*(-x2+x3))/(s1*s2)+((y1-y2)*(-y2+y3))/(s1*s2)+((z1-z2)*(-z2+z3))/(s1*s2)))*((x1-x2)/(s1*s2)-((x1-x2)*(-x2+x3)*std::abs(-x2+x3)*AbsDer(-x2+x3))/(s1*s12)-((y1-y2)*(-y2+y3)*std::abs(-x2+x3)*AbsDer(-x2+x3))/(s1*s12)-((z1-z2)*(-z2+z3)*std::abs(-x2+x3)*AbsDer(-x2+x3))/(s1*s12)))/(w2*std::sqrt(1-std::pow(((x1-x2)*(-x2+x3))/(s1*s2)+((y1-y2)*(-y2+y3))/(s1*s2)+((z1-z2)*(-z2+z3))/(s1*s2),2)));
	g(7)=(-2*(-b3+std::acos(((x1-x3)*(x2-x3))/(s4*s5)+((y1-y3)*(y2-y3))/(s4*s5)+((z1-z3)*(z2-z3))/(s4*s5)))*(-((y1-y3)/(s4*s5))-(y2-y3)/(s4*s5)+((x1-x3)*(x2-x3)*std::abs(y1-y3)*AbsDer(y1-y3))/(s6*s5)+((y1-y3)*(y2-y3)*std::abs(y1-y3)*AbsDer(y1-y3))/(s6*s5)+((z1-z3)*(z2-z3)*std::abs(y1-y3)*AbsDer(y1-y3))/(s6*s5)+((x1-x3)*(x2-x3)*std::abs(y2-y3)*AbsDer(y2-y3))/(s4*s10)+((y1-y3)*(y2-y3)*std::abs(y2-y3)*AbsDer(y2-y3))/(s4*s10)+((z1-z3)*(z2-z3)*std::abs(y2-y3)*AbsDer(y2-y3))/(s4*s10)))/(w3*std::sqrt(1-std::pow(((x1-x3)*(x2-x3))/(s4*s5)+((y1-y3)*(y2-y3))/(s4*s5)+((z1-z3)*(z2-z3))/(s4*s5),2)))-(2*(-b1+std::acos(((-x1+x2)*(-x1+x3))/(s7*s8)+((-y1+y2)*(-y1+y3))/(s7*s8)+((-z1+z2)*(-z1+z3))/(s7*s8)))*((-y1+y2)/(s7*s8)-((-x1+x2)*(-x1+x3)*std::abs(-y1+y3)*AbsDer(-y1+y3))/(s7*s11)-((-y1+y2)*(-y1+y3)*std::abs(-y1+y3)*AbsDer(-y1+y3))/(s7*s11)-((-z1+z2)*(-z1+z3)*std::abs(-y1+y3)*AbsDer(-y1+y3))/(s7*s11)))/(w1*std::sqrt(1-std::pow(((-x1+x2)*(-x1+x3))/(s7*s8)+((-y1+y2)*(-y1+y3))/(s7*s8)+((-z1+z2)*(-z1+z3))/(s7*s8),2)))-(2*(-b2+std::acos(((x1-x2)*(-x2+x3))/(s1*s2)+((y1-y2)*(-y2+y3))/(s1*s2)+((z1-z2)*(-z2+z3))/(s1*s2)))*((y1-y2)/(s1*s2)-((x1-x2)*(-x2+x3)*std::abs(-y2+y3)*AbsDer(-y2+y3))/(s1*s12)-((y1-y2)*(-y2+y3)*std::abs(-y2+y3)*AbsDer(-y2+y3))/(s1*s12)-((z1-z2)*(-z2+z3)*std::abs(-y2+y3)*AbsDer(-y2+y3))/(s1*s12)))/(w2*std::sqrt(1-std::pow(((x1-x2)*(-x2+x3))/(s1*s2)+((y1-y2)*(-y2+y3))/(s1*s2)+((z1-z2)*(-z2+z3))/(s1*s2),2)));
	g(8)=(-2*(-b3+std::acos(((x1-x3)*(x2-x3))/(s4*s5)+((y1-y3)*(y2-y3))/(s4*s5)+((z1-z3)*(z2-z3))/(s4*s5)))*(-((z1-z3)/(s4*s5))-(z2-z3)/(s4*s5)+((x1-x3)*(x2-x3)*std::abs(z1-z3)*AbsDer(z1-z3))/(s6*s5)+((y1-y3)*(y2-y3)*std::abs(z1-z3)*AbsDer(z1-z3))/(s6*s5)+((z1-z3)*(z2-z3)*std::abs(z1-z3)*AbsDer(z1-z3))/(s6*s5)+((x1-x3)*(x2-x3)*std::abs(z2-z3)*AbsDer(z2-z3))/(s4*s10)+((y1-y3)*(y2-y3)*std::abs(z2-z3)*AbsDer(z2-z3))/(s4*s10)+((z1-z3)*(z2-z3)*std::abs(z2-z3)*AbsDer(z2-z3))/(s4*s10)))/(w3*std::sqrt(1-std::pow(((x1-x3)*(x2-x3))/(s4*s5)+((y1-y3)*(y2-y3))/(s4*s5)+((z1-z3)*(z2-z3))/(s4*s5),2)))-(2*(-b1+std::acos(((-x1+x2)*(-x1+x3))/(s7*s8)+((-y1+y2)*(-y1+y3))/(s7*s8)+((-z1+z2)*(-z1+z3))/(s7*s8)))*((-z1+z2)/(s7*s8)-((-x1+x2)*(-x1+x3)*std::abs(-z1+z3)*AbsDer(-z1+z3))/(s7*s11)-((-y1+y2)*(-y1+y3)*std::abs(-z1+z3)*AbsDer(-z1+z3))/(s7*s11)-((-z1+z2)*(-z1+z3)*std::abs(-z1+z3)*AbsDer(-z1+z3))/(s7*s11)))/(w1*std::sqrt(1-std::pow(((-x1+x2)*(-x1+x3))/(s7*s8)+((-y1+y2)*(-y1+y3))/(s7*s8)+((-z1+z2)*(-z1+z3))/(s7*s8),2)))-(2*(-b2+std::acos(((x1-x2)*(-x2+x3))/(s1*s2)+((y1-y2)*(-y2+y3))/(s1*s2)+((z1-z2)*(-z2+z3))/(s1*s2)))*((z1-z2)/(s1*s2)-((x1-x2)*(-x2+x3)*std::abs(-z2+z3)*AbsDer(-z2+z3))/(s1*s12)-((y1-y2)*(-y2+y3)*std::abs(-z2+z3)*AbsDer(-z2+z3))/(s1*s12)-((z1-z2)*(-z2+z3)*std::abs(-z2+z3)*AbsDer(-z2+z3))/(s1*s12)))/(w2*std::sqrt(1-std::pow(((x1-x2)*(-x2+x3))/(s1*s2)+((y1-y2)*(-y2+y3))/(s1*s2)+((z1-z2)*(-z2+z3))/(s1*s2),2)));

	return g;
}

Matrix IsometricEnergyTriangleEvaluator::Hessian(const Attribute& attribute) const {
	return Matrix(0, 0);
}

static boost::tuples::tuple<Vec2n, Vec2n, Vec2n> TransformTo2D(const Vec3n& a, const Vec3n& b, const Vec3n& c, const Mat4x4n& T) {
	const Vec4n a_4 = T * math::V4n(a, 1.0);
	const Vec4n b_4 = T * math::V4n(b, 1.0);
	const Vec4n c_4 = T * math::V4n(c, 1.0);

	return boost::tuples::tuple<Vec2n, Vec2n, Vec2n>(math::V2n(a_4[0], a_4[1]), math::V2n(b_4[0], b_4[1]), math::V2n(c_4[0], c_4[1]));
}

nfloat DeformationEnergyTriangleEvaluator::f(const Attribute& attribute) const {
	const TriangleAttribute& ta = static_cast<const TriangleAttribute&>(attribute);
	STriangle* t = ta.t;

	const Mat4x4n T2 = t->Transform("position0");
	const Vec3n& a2 = t->vertex[0]->Get<Vec3n>("position0");
	const Vec3n& b2 = t->vertex[1]->Get<Vec3n>("position0");
	const Vec3n& c2 = t->vertex[2]->Get<Vec3n>("position0");
	boost::tuples::tuple<Vec2n, Vec2n, Vec2n> trans2 = TransformTo2D(a2, b2, c2, T2);
	const Vec2n P0 = boost::tuples::get<0>(trans2);
	const Vec2n P1 = boost::tuples::get<1>(trans2);
	const Vec2n P2 = boost::tuples::get<2>(trans2);

	Mat4x4n T1 = t->Transform(ta.strAttribute);
	const Vec3n& a1 = t->vertex[0]->Get<Vec3n>(ta.strAttribute);
	const Vec3n& b1 = t->vertex[1]->Get<Vec3n>(ta.strAttribute);
	const Vec3n& c1 = t->vertex[2]->Get<Vec3n>(ta.strAttribute);
	boost::tuples::tuple<Vec2n, Vec2n, Vec2n> trans1 = TransformTo2D(a1, b1, c1, T1);
	const Vec2n p0 = boost::tuples::get<0>(trans1);
	const Vec2n p1 = boost::tuples::get<1>(trans1);
	const Vec2n p2 = boost::tuples::get<2>(trans1);

	const nfloat& x0 = p0[0];
	const nfloat& x1 = p1[0];
	const nfloat& x2 = p2[0];
	const nfloat& y0 = p0[1];
	const nfloat& y1 = p1[1];
	const nfloat& y2 = p2[1];

	const nfloat& X0 = P0[0];
	const nfloat& X1 = P1[0];
	const nfloat& X2 = P2[0];
	const nfloat& Y0 = P0[1];
	const nfloat& Y1 = P1[1];
	const nfloat& Y2 = P2[1];

	const nfloat E = std::sqrt((std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2))*(std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)));
	return E - 2.0;
}

Vector DeformationEnergyTriangleEvaluator::Gradient(const Attribute& attribute) const {
	const TriangleAttribute& ta = static_cast<const TriangleAttribute&>(attribute);
	STriangle* t = ta.t;

	const Mat4x4n T2 = t->Transform("position0");
	const Vec3n& a2 = t->vertex[0]->Get<Vec3n>("position0");
	const Vec3n& b2 = t->vertex[1]->Get<Vec3n>("position0");
	const Vec3n& c2 = t->vertex[2]->Get<Vec3n>("position0");
	boost::tuples::tuple<Vec2n, Vec2n, Vec2n> trans2 = TransformTo2D(a2, b2, c2, T2);
	const Vec2n P0 = boost::tuples::get<0>(trans2);
	const Vec2n P1 = boost::tuples::get<1>(trans2);
	const Vec2n P2 = boost::tuples::get<2>(trans2);

	Mat4x4n T1 = t->Transform(ta.strAttribute);
	const Vec3n& a1 = t->vertex[0]->Get<Vec3n>(ta.strAttribute);
	const Vec3n& b1 = t->vertex[1]->Get<Vec3n>(ta.strAttribute);
	const Vec3n& c1 = t->vertex[2]->Get<Vec3n>(ta.strAttribute);
	boost::tuples::tuple<Vec2n, Vec2n, Vec2n> trans1 = TransformTo2D(a1, b1, c1, T1);
	const Vec2n p0 = boost::tuples::get<0>(trans1);
	const Vec2n p1 = boost::tuples::get<1>(trans1);
	const Vec2n p2 = boost::tuples::get<2>(trans1);
		
	const nfloat& x0 = p0[0];
	const nfloat& x1 = p1[0];
	const nfloat& x2 = p2[0];
	const nfloat& y0 = p0[1];
	const nfloat& y1 = p1[1];
	const nfloat& y2 = p2[1];

	const nfloat& X0 = P0[0];
	const nfloat& X1 = P1[0];
	const nfloat& X2 = P2[0];
	const nfloat& Y0 = P0[1];
	const nfloat& Y1 = P1[1];
	const nfloat& Y2 = P2[1];
	
	Vector g(9);
	g(0) = ((2*(std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2))*((X1-X2)*std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2))*AbsDer((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2))+(Y1-Y2)*std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2))*AbsDer((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2))))/(X2*(Y0-Y1)+X0*(Y1-Y2)+X1*(-Y0+Y2))+(2*(std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2))*((x1-x2)*(X2*(y0-y1)+X0*(y1-y2)+X1*(-y0+y2))*std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))*AbsDer((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))+(y1-y2)*(X2*(-y0+y1)+X1*(y0-y2)+X0*(-y1+y2))*std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))*AbsDer((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))+(x1-x2)*(Y0*(y1-y2)+Y1*y2-y1*Y2+y0*(-Y1+Y2))*std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))*AbsDer((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))+(y1-y2)*(Y0*(y1-y2)+Y1*y2-y1*Y2+y0*(-Y1+Y2))*std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))*AbsDer((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))))/std::pow(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2),2))/(2.*std::sqrt((std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2))*(std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2))));
	g(1) = ((2*(std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2))*(-((x1-x2)*(X0*(x1-x2)+X1*x2-x1*X2+x0*(-X1+X2))*std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))*AbsDer((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)))+(X0*(x1-x2)+X1*x2-x1*X2+x0*(-X1+X2))*(y1-y2)*std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))*AbsDer((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))-(x1-x2)*(x2*(-Y0+Y1)+x1*(Y0-Y2)+x0*(-Y1+Y2))*std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))*AbsDer((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))+(y1-y2)*(x2*(Y0-Y1)+x0*(Y1-Y2)+x1*(-Y0+Y2))*std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))*AbsDer((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))))/std::pow(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2),2)+(2*(std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2))*((X1-X2)*std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2))*AbsDer((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2))+(-Y1+Y2)*std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2))*AbsDer((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2))))/(X2*(-Y0+Y1)+X1*(Y0-Y2)+X0*(-Y1+Y2)))/(2.*std::sqrt((std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2))*(std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2))));
	g(2) = 0.0;
	g(3) = ((2*(std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2))*((-X0+X2)*std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2))*AbsDer((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2))+(-Y0+Y2)*std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2))*AbsDer((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2))))/(X2*(Y0-Y1)+X0*(Y1-Y2)+X1*(-Y0+Y2))+(2*(std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2))*((x0-x2)*(X2*(-y0+y1)+X1*(y0-y2)+X0*(-y1+y2))*std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))*AbsDer((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))+(y0-y2)*(X2*(y0-y1)+X0*(y1-y2)+X1*(-y0+y2))*std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))*AbsDer((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))+(x0-x2)*(-(Y1*y2)+Y0*(-y1+y2)+y0*(Y1-Y2)+y1*Y2)*std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))*AbsDer((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))+(y0-y2)*(-(Y1*y2)+Y0*(-y1+y2)+y0*(Y1-Y2)+y1*Y2)*std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))*AbsDer((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))))/std::pow(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2),2))/(2.*std::sqrt((std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2))*(std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2))));
	g(4) = ((2*(std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2))*(-((x0-x2)*(-(X1*x2)+X0*(-x1+x2)+x0*(X1-X2)+x1*X2)*std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))*AbsDer((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)))-(X0*(x1-x2)+X1*x2-x1*X2+x0*(-X1+X2))*(y0-y2)*std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))*AbsDer((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))-(x0-x2)*(x2*(Y0-Y1)+x0*(Y1-Y2)+x1*(-Y0+Y2))*std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))*AbsDer((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))+(y0-y2)*(x2*(-Y0+Y1)+x1*(Y0-Y2)+x0*(-Y1+Y2))*std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))*AbsDer((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))))/std::pow(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2),2)+(2*(std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2))*((X0-X2)*std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2))*AbsDer((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2))+(-Y0+Y2)*std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2))*AbsDer((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2))))/(X2*(Y0-Y1)+X0*(Y1-Y2)+X1*(-Y0+Y2)))/(2.*std::sqrt((std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2))*(std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2))));
	g(5) = 0.0;
	g(6) = ((2*(std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2))*((X0-X1)*std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2))*AbsDer((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2))+(Y0-Y1)*std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2))*AbsDer((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2))))/(X2*(Y0-Y1)+X0*(Y1-Y2)+X1*(-Y0+Y2))+(2*(std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2))*((x0-x1)*(X2*(y0-y1)+X0*(y1-y2)+X1*(-y0+y2))*std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))*AbsDer((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))+(y0-y1)*(X2*(-y0+y1)+X1*(y0-y2)+X0*(-y1+y2))*std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))*AbsDer((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))+(x0-x1)*(Y0*(y1-y2)+Y1*y2-y1*Y2+y0*(-Y1+Y2))*std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))*AbsDer((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))-(y0-y1)*(-(Y1*y2)+Y0*(-y1+y2)+y0*(Y1-Y2)+y1*Y2)*std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))*AbsDer((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))))/std::pow(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2),2))/(2.*std::sqrt((std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2))*(std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2))));
	g(7) = ((2*(std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2))*((x0-x1)*(-(X1*x2)+X0*(-x1+x2)+x0*(X1-X2)+x1*X2)*std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))*AbsDer((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))+(X0*(x1-x2)+X1*x2-x1*X2+x0*(-X1+X2))*(y0-y1)*std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))*AbsDer((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))+(x0-x1)*(x2*(Y0-Y1)+x0*(Y1-Y2)+x1*(-Y0+Y2))*std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))*AbsDer((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))+(y0-y1)*(x2*(Y0-Y1)+x0*(Y1-Y2)+x1*(-Y0+Y2))*std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))*AbsDer((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2))))/std::pow(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2),2)+(2*(std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2))*((-X0+X1)*std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2))*AbsDer((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2))+(Y0-Y1)*std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2))*AbsDer((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2))))/(X2*(Y0-Y1)+X0*(Y1-Y2)+X1*(-Y0+Y2)))/(2.*std::sqrt((std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(x1*y0-x2*y0-x0*y1+x2*y1+x0*y2-x1*y2)),2))*(std::pow(std::abs((X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((X1*y0-X2*y0-X0*y1+X2*y1+X0*y2-X1*y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2)+std::pow(std::abs((Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2)/(X1*Y0-X2*Y0-X0*Y1+X2*Y1+X0*Y2-X1*Y2)),2))));
	g(8) = 0.0;

	T1.Inverse();
	for(size_t i = 0; i < 3; i++) {
		Vec4n giT = T1 * math::V4n(g(3 * i), g(3 * i + 1), g(3 * i + 2), 0.0);
		g(3 * i) = giT[0];
		g(3 * i + 1) = giT[1];
		g(3 * i + 2) = giT[2];
	}

	return g;
}

Matrix DeformationEnergyTriangleEvaluator::Hessian(const Attribute& attribute) const {
	return Matrix(0, 0);
}

nfloat StretchEnergyTriangleEvaluator::f(const Attribute& attribute) const {
	const TriangleAttribute& ta = static_cast<const TriangleAttribute&>(attribute);
	STriangle* t = ta.t;

	const Vec3n& P0 = t->vertex[0]->Get<Vec3n>("position0");
	const Vec3n& P1 = t->vertex[1]->Get<Vec3n>("position0");
	const Vec3n& P2 = t->vertex[2]->Get<Vec3n>("position0");

	const Vec3n& p0 = t->vertex[0]->Get<Vec3n>(ta.strAttribute);
	const Vec3n& p1 = t->vertex[1]->Get<Vec3n>(ta.strAttribute);
	const Vec3n& p2 = t->vertex[2]->Get<Vec3n>(ta.strAttribute);

	const nfloat& x0 = p0[0];
	const nfloat& x1 = p1[0];
	const nfloat& x2 = p2[0];
	const nfloat& y0 = p0[1];
	const nfloat& y1 = p1[1];
	const nfloat& y2 = p2[1];

	const nfloat& X0 = P0[0];
	const nfloat& X1 = P1[0];
	const nfloat& X2 = P2[0];
	const nfloat& Y0 = P0[1];
	const nfloat& Y1 = P1[1];
	const nfloat& Y2 = P2[1];
	const nfloat& Z0 = P0[2];
	const nfloat& Z1 = P1[2];
	const nfloat& Z2 = P2[2];

	const nfloat E = std::sqrt((std::pow(X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2,2)+std::pow(X2*(-y0+y1)+X1*(y0-y2)+X0*(-y1+y2),2)+std::pow(x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2,2)+std::pow(Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2,2)+std::pow(x1*Z0-x2*Z0-x0*Z1+x2*Z1+x0*Z2-x1*Z2,2)+std::pow(y1*Z0-y2*Z0-y0*Z1+y2*Z1+y0*Z2-y1*Z2,2))/std::pow(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2),2))/Sqrt(2);

	return E;
}

Vector StretchEnergyTriangleEvaluator::Gradient(const Attribute& attribute) const {
	const TriangleAttribute& ta = static_cast<const TriangleAttribute&>(attribute);
	STriangle* t = ta.t;

	const Vec3n& P0 = t->vertex[0]->Get<Vec3n>("position0");
	const Vec3n& P1 = t->vertex[1]->Get<Vec3n>("position0");
	const Vec3n& P2 = t->vertex[2]->Get<Vec3n>("position0");

	const Vec3n& p0 = t->vertex[0]->Get<Vec3n>(ta.strAttribute);
	const Vec3n& p1 = t->vertex[1]->Get<Vec3n>(ta.strAttribute);
	const Vec3n& p2 = t->vertex[2]->Get<Vec3n>(ta.strAttribute);

	const nfloat& x0 = p0[0];
	const nfloat& x1 = p1[0];
	const nfloat& x2 = p2[0];
	const nfloat& y0 = p0[1];
	const nfloat& y1 = p1[1];
	const nfloat& y2 = p2[1];

	const nfloat& X0 = P0[0];
	const nfloat& X1 = P1[0];
	const nfloat& X2 = P2[0];
	const nfloat& Y0 = P0[1];
	const nfloat& Y1 = P1[1];
	const nfloat& Y2 = P2[1];
	const nfloat& Z0 = P0[2];
	const nfloat& Z1 = P1[2];
	const nfloat& Z2 = P2[2];
	
	const nfloat s0 = std::pow(X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2,2);
	const nfloat s1 = std::pow(x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2,2);
	const nfloat s2 = std::pow(x1*Z0-x2*Z0-x0*Z1+x2*Z1+x0*Z2-x1*Z2,2);
	const nfloat s3 = std::pow(X2*(-y0+y1)+X1*(y0-y2)+X0*(-y1+y2),2);
	const nfloat s4 = std::pow(Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2,2);
	const nfloat s5 = std::pow(y1*Z0-y2*Z0-y0*Z1+y2*Z1+y0*Z2-y1*Z2,2);
	const nfloat s6 = std::pow(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2),3);
	const nfloat s7 = std::pow(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2),2);
	const nfloat s8 = std::sqrt(2.0);
	const nfloat s9 = std::sqrt((s0+s3+s1+s4+s2+s5)/s7);
	const nfloat s10 = 1.0/(2.*s8*s6*s9);
	
	Vector g(9);
	g(0) = (-2*(-y1+y2)*(s0+s1+s2)-2*(-y1+y2)*(s3+s4+s5)+2*(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2))*((X1-X2)*(-(X1*x2)+X0*(-x1+x2)+x0*(X1-X2)+x1*X2)+(-Y1+Y2)*(x2*(-Y0+Y1)+x1*(Y0-Y2)+x0*(-Y1+Y2))+(-Z1+Z2)*(x2*(-Z0+Z1)+x1*(Z0-Z2)+x0*(-Z1+Z2))))*s10;
	g(1) = (-2*(x1-x2)*(s0+s1+s2)-2*(x1-x2)*(s3+s4+s5)+2*(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2))*((X1-X2)*(X2*(-y0+y1)+X1*(y0-y2)+X0*(-y1+y2))+(Y1-Y2)*(-(Y1*y2)+Y0*(-y1+y2)+y0*(Y1-Y2)+y1*Y2)+(-Z1+Z2)*(y2*(-Z0+Z1)+y1*(Z0-Z2)+y0*(-Z1+Z2))))*s10;
	g(2) = 0.0;
	g(3) = (-2*(y0-y2)*(s0+s1+s2)-2*(y0-y2)*(s3+s4+s5)+2*(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2))*((X0-X2)*(X0*(x1-x2)+X1*x2-x1*X2+x0*(-X1+X2))+(Y0-Y2)*(x2*(-Y0+Y1)+x1*(Y0-Y2)+x0*(-Y1+Y2))+(Z0-Z2)*(x2*(-Z0+Z1)+x1*(Z0-Z2)+x0*(-Z1+Z2))))*s10;
	g(4) = (-2*(-x0+x2)*(s0+s1+s2)-2*(-x0+x2)*(s3+s4+s5)+2*(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2))*((X0-X2)*(X2*(y0-y1)+X0*(y1-y2)+X1*(-y0+y2))+(Y0-Y2)*(Y0*(y1-y2)+Y1*y2-y1*Y2+y0*(-Y1+Y2))+(Z0-Z2)*(y2*(-Z0+Z1)+y1*(Z0-Z2)+y0*(-Z1+Z2))))*s10;
	g(5) = 0.0;
	g(6) = (-2*(-y0+y1)*(s0+s1+s2)-2*(-y0+y1)*(s3+s4+s5)+2*(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2))*((-X0+X1)*(X0*(x1-x2)+X1*x2-x1*X2+x0*(-X1+X2))+(Y0-Y1)*(x2*(Y0-Y1)+x0*(Y1-Y2)+x1*(-Y0+Y2))+(Z0-Z1)*(x2*(Z0-Z1)+x0*(Z1-Z2)+x1*(-Z0+Z2))))*s10;
	g(7) = (-2*(x0-x1)*(s0+s1+s2)-2*(x0-x1)*(s3+s4+s5)+2*(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2))*((X0-X1)*(X2*(-y0+y1)+X1*(y0-y2)+X0*(-y1+y2))+(-Y0+Y1)*(Y0*(y1-y2)+Y1*y2-y1*Y2+y0*(-Y1+Y2))+(Z0-Z1)*(y2*(Z0-Z1)+y0*(Z1-Z2)+y1*(-Z0+Z2))))*s10;
	g(8) = 0.0;
	
// 	Vector g(9);
// 	g(0) = (-2*(-y1+y2)*(std::pow(X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2,2)+std::pow(x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2,2)+std::pow(x1*Z0-x2*Z0-x0*Z1+x2*Z1+x0*Z2-x1*Z2,2))-2*(-y1+y2)*(std::pow(X2*(-y0+y1)+X1*(y0-y2)+X0*(-y1+y2),2)+std::pow(Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2,2)+std::pow(y1*Z0-y2*Z0-y0*Z1+y2*Z1+y0*Z2-y1*Z2,2))+2*(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2))*((X1-X2)*(-(X1*x2)+X0*(-x1+x2)+x0*(X1-X2)+x1*X2)+(-Y1+Y2)*(x2*(-Y0+Y1)+x1*(Y0-Y2)+x0*(-Y1+Y2))+(-Z1+Z2)*(x2*(-Z0+Z1)+x1*(Z0-Z2)+x0*(-Z1+Z2))))/(2.*Sqrt(2)*std::pow(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2),3)*Sqrt((std::pow(X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2,2)+std::pow(X2*(-y0+y1)+X1*(y0-y2)+X0*(-y1+y2),2)+std::pow(x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2,2)+std::pow(Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2,2)+std::pow(x1*Z0-x2*Z0-x0*Z1+x2*Z1+x0*Z2-x1*Z2,2)+std::pow(y1*Z0-y2*Z0-y0*Z1+y2*Z1+y0*Z2-y1*Z2,2))/std::pow(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2),2)));
// 	g(1) = (-2*(x1-x2)*(std::pow(X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2,2)+std::pow(x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2,2)+std::pow(x1*Z0-x2*Z0-x0*Z1+x2*Z1+x0*Z2-x1*Z2,2))-2*(x1-x2)*(std::pow(X2*(-y0+y1)+X1*(y0-y2)+X0*(-y1+y2),2)+std::pow(Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2,2)+std::pow(y1*Z0-y2*Z0-y0*Z1+y2*Z1+y0*Z2-y1*Z2,2))+2*(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2))*((X1-X2)*(X2*(-y0+y1)+X1*(y0-y2)+X0*(-y1+y2))+(Y1-Y2)*(-(Y1*y2)+Y0*(-y1+y2)+y0*(Y1-Y2)+y1*Y2)+(-Z1+Z2)*(y2*(-Z0+Z1)+y1*(Z0-Z2)+y0*(-Z1+Z2))))/(2.*Sqrt(2)*std::pow(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2),3)*Sqrt((std::pow(X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2,2)+std::pow(X2*(-y0+y1)+X1*(y0-y2)+X0*(-y1+y2),2)+std::pow(x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2,2)+std::pow(Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2,2)+std::pow(x1*Z0-x2*Z0-x0*Z1+x2*Z1+x0*Z2-x1*Z2,2)+std::pow(y1*Z0-y2*Z0-y0*Z1+y2*Z1+y0*Z2-y1*Z2,2))/std::pow(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2),2)));
// 	g(2) = 0.0;
// 	g(3) = (-2*(y0-y2)*(std::pow(X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2,2)+std::pow(x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2,2)+std::pow(x1*Z0-x2*Z0-x0*Z1+x2*Z1+x0*Z2-x1*Z2,2))-2*(y0-y2)*(std::pow(X2*(-y0+y1)+X1*(y0-y2)+X0*(-y1+y2),2)+std::pow(Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2,2)+std::pow(y1*Z0-y2*Z0-y0*Z1+y2*Z1+y0*Z2-y1*Z2,2))+2*(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2))*((X0-X2)*(X0*(x1-x2)+X1*x2-x1*X2+x0*(-X1+X2))+(Y0-Y2)*(x2*(-Y0+Y1)+x1*(Y0-Y2)+x0*(-Y1+Y2))+(Z0-Z2)*(x2*(-Z0+Z1)+x1*(Z0-Z2)+x0*(-Z1+Z2))))/(2.*Sqrt(2)*std::pow(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2),3)*Sqrt((std::pow(X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2,2)+std::pow(X2*(-y0+y1)+X1*(y0-y2)+X0*(-y1+y2),2)+std::pow(x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2,2)+std::pow(Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2,2)+std::pow(x1*Z0-x2*Z0-x0*Z1+x2*Z1+x0*Z2-x1*Z2,2)+std::pow(y1*Z0-y2*Z0-y0*Z1+y2*Z1+y0*Z2-y1*Z2,2))/std::pow(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2),2)));
// 	g(4) = (-2*(-x0+x2)*(std::pow(X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2,2)+std::pow(x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2,2)+std::pow(x1*Z0-x2*Z0-x0*Z1+x2*Z1+x0*Z2-x1*Z2,2))-2*(-x0+x2)*(std::pow(X2*(-y0+y1)+X1*(y0-y2)+X0*(-y1+y2),2)+std::pow(Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2,2)+std::pow(y1*Z0-y2*Z0-y0*Z1+y2*Z1+y0*Z2-y1*Z2,2))+2*(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2))*((X0-X2)*(X2*(y0-y1)+X0*(y1-y2)+X1*(-y0+y2))+(Y0-Y2)*(Y0*(y1-y2)+Y1*y2-y1*Y2+y0*(-Y1+Y2))+(Z0-Z2)*(y2*(-Z0+Z1)+y1*(Z0-Z2)+y0*(-Z1+Z2))))/(2.*Sqrt(2)*std::pow(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2),3)*Sqrt((std::pow(X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2,2)+std::pow(X2*(-y0+y1)+X1*(y0-y2)+X0*(-y1+y2),2)+std::pow(x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2,2)+std::pow(Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2,2)+std::pow(x1*Z0-x2*Z0-x0*Z1+x2*Z1+x0*Z2-x1*Z2,2)+std::pow(y1*Z0-y2*Z0-y0*Z1+y2*Z1+y0*Z2-y1*Z2,2))/std::pow(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2),2)));
// 	g(5) = 0.0;
// 	g(6) = (-2*(-y0+y1)*(std::pow(X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2,2)+std::pow(x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2,2)+std::pow(x1*Z0-x2*Z0-x0*Z1+x2*Z1+x0*Z2-x1*Z2,2))-2*(-y0+y1)*(std::pow(X2*(-y0+y1)+X1*(y0-y2)+X0*(-y1+y2),2)+std::pow(Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2,2)+std::pow(y1*Z0-y2*Z0-y0*Z1+y2*Z1+y0*Z2-y1*Z2,2))+2*(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2))*((-X0+X1)*(X0*(x1-x2)+X1*x2-x1*X2+x0*(-X1+X2))+(Y0-Y1)*(x2*(Y0-Y1)+x0*(Y1-Y2)+x1*(-Y0+Y2))+(Z0-Z1)*(x2*(Z0-Z1)+x0*(Z1-Z2)+x1*(-Z0+Z2))))/(2.*Sqrt(2)*std::pow(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2),3)*Sqrt((std::pow(X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2,2)+std::pow(X2*(-y0+y1)+X1*(y0-y2)+X0*(-y1+y2),2)+std::pow(x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2,2)+std::pow(Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2,2)+std::pow(x1*Z0-x2*Z0-x0*Z1+x2*Z1+x0*Z2-x1*Z2,2)+std::pow(y1*Z0-y2*Z0-y0*Z1+y2*Z1+y0*Z2-y1*Z2,2))/std::pow(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2),2)));
// 	g(7) = (-2*(x0-x1)*(std::pow(X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2,2)+std::pow(x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2,2)+std::pow(x1*Z0-x2*Z0-x0*Z1+x2*Z1+x0*Z2-x1*Z2,2))-2*(x0-x1)*(std::pow(X2*(-y0+y1)+X1*(y0-y2)+X0*(-y1+y2),2)+std::pow(Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2,2)+std::pow(y1*Z0-y2*Z0-y0*Z1+y2*Z1+y0*Z2-y1*Z2,2))+2*(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2))*((X0-X1)*(X2*(-y0+y1)+X1*(y0-y2)+X0*(-y1+y2))+(-Y0+Y1)*(Y0*(y1-y2)+Y1*y2-y1*Y2+y0*(-Y1+Y2))+(Z0-Z1)*(y2*(Z0-Z1)+y0*(Z1-Z2)+y1*(-Z0+Z2))))/(2.*Sqrt(2)*std::pow(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2),3)*Sqrt((std::pow(X0*x1-x0*X1-X0*x2+X1*x2+x0*X2-x1*X2,2)+std::pow(X2*(-y0+y1)+X1*(y0-y2)+X0*(-y1+y2),2)+std::pow(x1*Y0-x2*Y0-x0*Y1+x2*Y1+x0*Y2-x1*Y2,2)+std::pow(Y0*y1-y0*Y1-Y0*y2+Y1*y2+y0*Y2-y1*Y2,2)+std::pow(x1*Z0-x2*Z0-x0*Z1+x2*Z1+x0*Z2-x1*Z2,2)+std::pow(y1*Z0-y2*Z0-y0*Z1+y2*Z1+y0*Z2-y1*Z2,2))/std::pow(x2*(-y0+y1)+x1*(y0-y2)+x0*(-y1+y2),2)));
// 	g(8) = 0.0;

	return g;
}

Matrix StretchEnergyTriangleEvaluator::Hessian(const Attribute& attribute) const {
	return Matrix(0, 0);
}

}
